Draft version February 2, 2008 

Preprint typeset using LATp^X style emulateapj v. 04/03/99 



SIMULATING WEAK LENSING BY LARGE SCALE STRUCTURE 



O 
O 



in 



> 
in 
in 
in 
m 
O 

o 

Or 
i 

o 

c3 



13 



Chris Vale and Martin White 

Departments of Astronomy and Physics, University of California, Berkeley, CA 94720 
Draft version February 2, 2008 

ABSTRACT 

We model weak gravitational lensing of light by large-scale structure using ray tracing through N-body 
simulations. The method is described with particular attention paid to numerical convergence. We 
investigate some of the key approximations in the multi-plane ray tracing algorithm. Our simulated 
shear and convergence maps are used to explore how well standard assumptions about weak lensing 
hold, especially near large peaks in the lensing signal. 

Subject headings: weak lensing, cosmology-theory 



1. INTRODUCTION 

Weak gravitational lensing is becoming an indispens- 
able tool in modern cosmology. Building on early work 
by Gunnjgp7|), Blandford et al. l(TWT|> and Miralda- 
Escude (1991) showed that the shear and magnification are 
on the order of a few percent in popular cosmologies, which 
is now within reach of observations. Indeed detections of 
shear correlations by several groups (Bacon et al. [2000 
Van Waerbeke et al. I2000al R hodes et al. I2UUT1 Hoekstra 
et al. 11002 Brown et al.EOOl Jarvis et al. l!-003"|) indicate 
the dramatic progress being made in observational power 
and precision. Other projects, such as the CFHT legacy 
survey 1 , the Deep Lens Survey 2 , and NO AO deep survey 3 
and proposed projects such as SNAP 4 , DMT/LSST 5 , and 
Pan— STARRS 6 will continue this trend as they map the 
shear on large fractions of the sky. In order for the full 
scientific return from lensing to be realized, both funda- 
mental theory and data analysis techniques must keep pace 
with these observational advances. 

In principle the effect exploited by these surveys is sim- 
ple: gravitational lensing of light rays by foreground large- 
scale structure magnifies and shears the images of back- 
ground galaxies, and this distortion can be mapped to pro- 
vide information on the matter distribution and cosmologi- 
cal model. Like the cosmic microwave background (CMB), 
the theory of weak lensing is simple and clean (see Bartel- 
mann & Schneider 2001 or Mellier 1999 for a review). Un- 
like the the CMB, however, lensing on sub-degree scales 
probes large-scale structure in the non-linear regime and a 
full description has not been achieved through purely an- 
alytic means. Fortunately, on scales relevant to weak lens- 
ing, the hierarchical growth of structure can be accurately 
simulated using N-body methods. The lensing effect of 
this structure can then be computed directly by ray tracing 
through the simulation. In this paper we describe an im- 
plementation of such a method, building upon the work of 
(Blandford et aL llMTl Wambsg anss, C en & Ostriker lT^gl 
Couchman, Barber, & Thomas 119981 Fluke, Webster, & 

1 http: / / www.cfht . hawaii . edu / Science /CFHLS 
2 http://dls. bell-labs, com/ 
3 http: / / www.noao.edu/noao / noaodeep / 
4 http:/ /snap. lbl.gov 

5 http: / /www. dmtelescope.org/dark_home. html 
6 http: // www. ifa. hawaii . edu / pan-starrs / 



Mortlock ITT)9"gl Hama na, Martel fc Futa mase 1211(101 Jain, 
Seljak, & White l2li001 Whi te & Hu 2000!; Barb er, Th omas, 
Couchman, fc Flu ke 1217001 Hamana & Mellier E_00T1 Hen - 
nawi et al. 1217011 Padmanabhan, Seljak, & Pen l!_00!_|) . 

Our goal is to quantify the numerical precision of such 
algorithms, and to investigate the computational cost of 
achieving a given theoretical precision. 

The outline of the paper is as follows. In [|2]we introduce 
the weak lensing formalism and in iQ we describe our im- 
plementation of a multi-plane ray-tracing algorithm. In 21 
we present some basic results for comparison to previous 
work while in |S] we evaluate the numerical convergence of 
the algorithm. Some interesting physical results are given 
in S|0J before we conclude in SJ7| 

2. WEAK LENSING FORMALISM 

We make use of a multiple lens plane algorithm in order 
to simulate the distortion and magnification effect of fore- 
ground matter on background light rays. In the follow- 
ing section, we provide a summary of equations directly 
relevant to this paper (see Jain, Seljak, & White (2000) 
or Schneider, Ehlers, & Falco (1992) for a more thor- 
ough treatment) and then describe our ray tracing method. 
Note that we adopt units where c = 1, and we will work 
in comoving coordinates. 

In a universe governed by the Robertson- Walker metric, 
the change in direction of a light ray propagating through 
space is: 



da = -2V_l^ dx 



(1) 



where da is the bend angle, denotes the spatial gradi- 
ent perpendicular to the path of the light ray, <p is the 3 
dimensional peculiar gravitational potential, and x is the 
radial comoving coordinate. 

The change in position on a plane perpendicular to the 
line of sight at a position \ due to a deflection da at %' 
is: 



dx(x) = r(x - x')da{x) 



(2) 



where r(x) is the comoving angular diameter distance. In- 
tegrating to get x{x) and then dividing by r{x) yields: 



0(X) 



r(x) 



dx' r(x-x')Vx<t> + 6(0) 



(3) 
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where 6(x) i s the angular position of the light ray at 
comoving coordinate x- The 2x2 distortion matrix, 
Ay =^i(x)/9^(0), is given by: 



Aij — 2 



Jo 



<j) + Si. 



where we make use of the definition: 

r(x ~ X')r(x') 



g{x',x) 



ix) 



This matrix is normally decomposed as: 

A _ / 1 - k - 71 -72 - LJ 
-72 + W 1 — K + 71 



(4) 



(5) 



(G) 



where k is the convergence, 7 is the shear, and to is the 
rotation. 

The gravitational potential, </>', is related to the mass 
density, p, through Poisson's equation: 



72 jj 



4irGa 2 p 



(7) 



where the gradient is with respect to comoving coordinates 
(that's why there's an a 2 on the right hand side) and </>' is 
the sum of a smooth background potential <f> and the lo- 
cal peculiar potential (f>. We can subtract the background 
terms and use V 2 </> = AirGp, where p is the mean density 
in the region, and substitute a 3 p = po to find 



AnGpo ( p 



(8) 



If we then define the relative mass overdensity, 6 = p/p—1, 
and substitute p~o = 3HQfl m /8irG we obtain 



; 2 1 3-ff 2 f2 m <5 

~~ 2a 



(9) 



In the discrete approximation, the interval between the 
source and the observer is divided into distinct regions 
separated by a comoving distance A%. The relative mass 
overdensity in a given region is then projected onto a lens 
plane within the region and perpendicular to the line of 
sight. An effective two dimensional gravitational potential, 
ip = 2 J 4>dx, is defined on the this plane, enabling a two 
dimensional Poisson equation relating ip and the matter 
content in the region to be written as 



VV = 3H 2 n rn T, 



(10) 



where S = J S dx is the projected two dimensional relative 
mass overdensity. Eqs © and J3J then become 



r(x„ - Xm) 



A„ = I 



n-1 

E 

ni—l 



(ii) 



where U m is the shear tensor in the m th region defined by 

d 2 ip m 



dxidxj 



(12) 



It will be useful to decompose the distortion matrix at a 
given plane into two components (Seitz, Schneider, and 
Ehlcrs 1994) B„ and C„ such that for a zero curvature 
universe (assumed hereafter): 



A ra +i = I — B„ + C n 
B„ = B„_i + x n U„A„ 

v " -(C„_l + XnU„A n ) 



Xn+l 



3. RAY-TRACING ALGORITHM 



(13) 



The goal of our weak lensing algorithm is to simulate the 
deflection and distortion that light rays would experience 
as they propagate through an intervening matter distribu- 
tion that is statistically similar to that of the real universe. 
In this section we describe two different implementations 
of the weak lensing algorithm: first, a standard implemen- 
tation based on multiple projected-mass lens planes, and 
then a second, three-dimensional version. The modeling 
of the mass distribution is discussed in fj5] 

3.1. Standard (2D) algorithm 

The standard implementation consists of creating a two- 
dimensional projected mass plane from an N-body simu- 
lation, computing the position, propagation direction, and 
distortion matrix on the plane for each ray, and then re- 
peating the process until the rays have been traced from 
the observer to the source. Here is a more detailed de- 
scription: 

1. We begin by using an N-body simulation (which we 
describe later) to create the structure responsible for the 
lensing. Since lensing is only sensitive to density fluctu- 
ations perpendicular to the line of sight, the simulation 
must resolve structures on length scales ~ X 9, where x is 
a typical survey depth, and 9 is the desired angular reso- 
lution. Light-rays from rcdshifts of interest must therefore 
travel thousands of Mpc through structures resolved below 
Mpc scales. 

It is currently computationally very expensive to simu- 
late both scales simultaneously, but, fortunately, this is not 
necessary. Since there is little lensing power on very large 
scales, it makes no sense to simulate huge boxes whose 
pieces are almost independent. We employ a technique 
(common in the literature) in which we simulate a volume 
which is large enough to be a fair sample of the universe, 
yet small enough to provide sufficient small scale resolu- 
tion at reasonable cost. In this simulation, mass particles 
in a box evolve in time under the influence of gravity, and 
the positions of the particles are recorded at time intervals 
corresponding to regular intervals in comoving distance. 
These output boxes are then used to create the mass dis- 
tribution to be traced through. Note that since each box 
is actually the same matter distribution at different stages 
of evolution, care must be taken to avoid tracing over the 
same structures more than once. 

2. A selected number of light rays to be traced are then 
initialized at the observer for a square field of view. For 
each ray, we track the position x, the propagation direction 
a n relative to the line of sight, and the matrices B and C 
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from equation (|13|) . The recursion relations (Eq. Illjl for 
A„ and 9 n require as inputs the previous values for V 
and U, all of which need to be stored. However, we make 
use of a less memory intensive algorithm for the distortion 
matrix based on Eq. i|13|) in which we track the values the 
decomposed distortion matrix B and C for each ray. Each 
iteration requires as a new input only the current value 
of U at each ray location, which may then be discarded. 
Similarly, by tracking d? n , we require only one value of 
V^V a t a time. 

3. Next, we select the source redshift and the number of 
lens planes we wish to use. Note that the sources are all 
initialized at the same redshift, and the planes are equally 
spaced in comoving distance. This sets the comoving dis- 
tance between planes, A%. We create a lens plane from 
an N-body simulation box using the following technique. 
We select the N-body particle dump closest in comoving 
distance to that of the desired lens plane, and then ran- 
domly choose the x, y, or z axis as the propagation axis, 
X- A two-dimensional N x N lens plane grid is then created 
by projecting mass particles from the N-body simulation 
in a direction parallel to the \ ax i s onto a plane at the 
desired comoving distance, and then recording the mass 
onto the grid locations using a cloud in cell (CIC) assign- 
ment scheme. All of the particles in the box within 
of the desired plane are used to make the grid, which is 
then normalized by the mean density. The projected field 
is thus a square of side length equal to that of the box. 

The origin of the lens plane is randomly selected. This 
origin shift is facilitated by the fact that the N-body sim- 
ulation uses periodic boundary conditions, and we can use 
this periodicity to re-map the mass density about any ori- 
gin we choose. This process of rotation and origin shift 
makes it unlikely that the same structures will be traced 
repeatedly. Note that both the origin and orientation of 
the axes are maintained until the box is traced all the 
way through, at which point the box is re-randomized, so 
that consecutive lens planes are made from a matter dis- 
tribution that is continuous everywhere except at the box 
boundaries. 

A variation of this technique, to our knowledge pre- 
sented here for the first time, and which we make use 
of for comparison, avoids introducing these discontinuities 
entirely. The propagation direction is selected at an angle 
relative to one of the box axes. The rays are initialized 
about a central line of sight along this direction, and the 
lens planes are created perpendicular to the central line of 
sight and with their origins on it. The angle is chosen so 
that the rays, traced from the observer to the source, never 
propagate through the same volume twice. However, since 
the boxes are never rotated, nor is there any random origin 
shift, the lens planes are made from a matter distribution 
that is continuous everywhere. 

4. The two component gradient of the potential and 
the four components of the matrix U are then computed 
on the grid. This is done by first taking the Fourier trans- 
formation of the mass density plane, then taking the first 
or second derivatives, respectively, and computing the in- 
verse transformation. We do this by making use of the 
relation of real space derivatives in Fourier space, for ex- 



ample 

^ = -i [ k x me~^ s d 2 k (14) 

5. The position x n for each ray on the plane is computed 
using the stored values x n -\ and <3„_i. Values for U and 
Vi/> are then assigned at each ray position, again using 
a CIC assignment scheme, and the new deflection angle 
a„ = a„_i + Aa„ is then computed where 

Aa = -V^ (15) 

follows from Eq. The distortion matrix A„ is tem- 
porarily created from stored values of B„_i and C„_i and 
used along with U„ to create the new values B„ and C n 
for each ray on the plane. 

6. The process is repeated until each ray is traced all the 
way back to its source, thus ensuring that all the rays will 
meet at the observer, and the final angular position and 
distortion matrix are recorded for each ray. 

Clearly, one limitation of any algorithm based on lens 
planes is that the light rays themselves do not all prop- 
agate in precisely the same direction and cannot all be 
perpendicular to a given lens plane. The weak lensing for- 
malism of the previous section calls for the calculation of 
derivatives perpendicular to the light rays propagation di- 
rection, so the multiple lens plane algorithm is necessarily 
only approximate in this respect. 

3.2. Extended (3D) algorithm 

We have created a second version of the ray-tracing algo- 
rithm that is fully three-dimensional, and therefore doesn't 
suffer from this drawback. It is quite similar to the stan- 
dard implementation in most respects, and makes use of 
the same N-body particle dumps. Here is a brief descrip- 
tion. 

1. The source redshift is chosen and light rays in a square 
field of view are initialized at the observer. As before, we 
will track the position and propagation direction of each 
ray, and also the matrices B and C. Note that these are 
still 2x2 matrices, since, like the matrix A, they still 
describe the distortion of a two-dimensional image. 

2. The matter distribution from the particle dump nearest 
in redshift to the current position is mapped onto a three- 
dimensional regular grid using a CIC assignment scheme, 
and the grid is then normalized by the mean density. As 
before, it's important not to re-trace the same structures, 
so the box orientations and origins are randomly reas- 
signed, when necessary, to avoid this. 

3. The relationship in Fourier space between ip and S is 
used to calculate the three first derivatives and six in- 
dependent second derivatives of <fi at each point on the 
grid. These derivatives are then evaluated using the CIC 
method at each ray position. The components perpendicu- 
lar to each ray direction are calculated at each ray position 
and used to update the values of B, C and a n . The rays 
are then projected forward by a pre-selected comoving dis- 
tance, and the derivatives are then evaluated at the new 
positions. This process is repeated until it becomes time 
to use the next particle dump. 
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Fig. 1. — (left) An image of the convergence re from a single realization for a 2° field of view. The greyscale is linear in re, ranging from 
white (re = —0.04) to black (re > 0.7). Regions of high convergence are due to massive structures (typically clusters of galaxies) along the 
line-of-sight. (right) The shear field 7 for the same realization, exaggerated for clarity (the magnitude of the shear is at the percent level). 
Note that the shear is tangential around regions of high convergence. 



4. BASIC RESULTS 

Before we discuss the numerical convergence of our tech- 
nique, we present some basic science results. These are 
broadly in agreement with similar studies presented by 
(Jain Selj ak, & White l27)0Til Hamana, Martel, & Futa- 
mase|2E3ni White & HuEED!) as we elucidate later. Fig.ffl 
is an example of the convergence (k) and shear (7) maps 
that result from the ray-tracing procedure. The maps are 
2° on a side and contain 2048 2 lines of sight. The average 
magnitude of the convergence and shear is about 1.8%, 
close to the level predicted by Blandford et al. ( Tl991[l and 
Miralda-Escude ( KWfy . 

In the following section we describe these fields. We 
begin by presenting the 1-point distribution, and demon- 
strate that there are not likely to be any unlensed images 
at high redshift for this cosmology. We then present the 
power spectrum, which we find to be consistent with the 
Born and Limber approximations. We then discuss the 
skewness (S3) and the kurtosis (S4), which are given as a 
function of angular smoothing scale. We conclude with a 
description of the convergence peaks in our maps. 

The lowest order statistic for the maps is the distribu- 
tion of the shear amplitude or convergence. We show two 
typical examples in Fig. [2] in terms of the magnification, 
which in the weak lensing limit is simply /i = 1 + 2k. 
This is compared to the 'universal' form of Wang, Holz, 
& Munshi ((5002) in The average magnitude of the 
convergence over 10 realizations is 0.018, consistent with 
the expectation of a weak lensing effect on the order of a 
few percent. The mean and variance of the convergence 
are -0.002 and 7.5 x 10 -4 , respectively. The minimum k 




Fig. 2. — (top) A typical magnification PDF from our simulation 
(squares) compared with the prediction (solid line) using the method 
of Wang, Holz, & Munshi 120021 . The curves typically agree to 
within ~ 5% near the peaks and ~ 20% in the high fj, tails, (bottom) 
A second example where the model PDF is shifted on the horizontal 
axis with respect to that of the simulation. This occurs frequently 
in our maps. The shape of the measured PDF is still well fit by the 
prediction of the model (solid line), as can be seen when this curve 
is shifted in fi (dashed line). 
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for any of our realizations was —0.049, which is greater 
than the 'empty beam' value (an empty beam is a path 
from the source to the observer which is sufficiently void 
of matter that gravitational interactions may be ignored) 
of —0.064 for our cosmology, consistent with the fact that 
all the light rays in our simulations are lensed to some 
extent. 

The shear correlations, from which cosmological infer- 
ences are usually drawn, are shown in Fig. |31 in terms of 
angular power spectra. Figure|3]shows the power spectrum 
of (k) vs. the spectra of the E-modc of the shear (7^) and 
of the rotation (a>) for one of our simulated maps. The 
shear E-mode is defined to be the curl free component of 
the shear tensor, and is easily calculated in Fourier space 
using 

= TjT-172 ( 16 ) 



We present it here in the form of White & Hu ( 2000 ) 



x > K y 



where the tilde denotes the Fourier transform. The 'excess' 
power in k and je at t ~ 200 is just a fluctuation in this 
particular map. 

One of the most important predictions of the weak lens- 
ing approximation is that the shear matrix A can be well 
approximated using the 2nd derivatives of a scalar poten- 
tial, as would be the case if only a single lens plane were 
used to compute it. Then A would be symmetric, lu would 
be zero, and je would equal k. Thus, in the weak lensing 
regime, the spectra of k and je are predicted to be nearly 
equal to each other, and to be much larger than that of 
to. We agree with earlier work (e.g. JSW) in verifying this 
basic assumption of weak lensing (Fig. [3J . 

Figure 0] shows the angular power spectrum of k (at 
various levels angular resolution) compared to the semi- 
analytic result which makes both the Born and Lim- 
ber approximations. This result was first derived by 
Kaiser IjKaiser 1992|l and extended by Jain & Seljak (JT997 ). 



icr 4 - 




Fig. 3. — The angular power spectra of the convergence K, curl-free 
shear 7^, and rotation lo (multiplied by 10 4 ) from one of our simu- 
lated maps. As predicted for the weak lensing regime, the spectrum 
of lu is much less than the spectra of k and 7^, which are nearly 
equal. 



x'd x ' 



g(x',x) 

a(x') 



A t 2 n (fc, a ) (17) 



where 3 (x', x) is defined in Eq. ©, A 2 n (fc) = k 3 P(k)/(2ir 2 ) 
is the contribution to the mass variance per logarithmic 
interval in wavenumber and A 2 (£) = £ 2 Ce/(2ir) is the 
contribution to K 2 ms per logarithmic interval in angular 
wavenumber I. 



The agreement between the numerical and semi-analytic 
predictions is quite good, but not perfect. One reason 
for this could be that in evaluating Eq. I|17|l we used the 
method of Peacock & Dodds l|1996() to compute the non- 
linear power spectrum P{k). As can be seen in Figure [3] 
where we compare this prediction with the measured 3- 
dimensional matter power spectrum of our N-body simu- 
lation, the Peacock & Dodds IT996[) fitting formula is not 
exact. If we replace the A 2 n (fc) in Eq. (|17|l with the values 
actually measured in our simulation, the discrepancy dis- 
appears. We conclude therefore that the Born and Limber 
approximations are consistent with the power spectrum 
results from our simulations (we shall return to this issue 
in 

The maps in Fig. ^ are noticeably non-Gaussian, with 
identifiable objects accounting for large positive values of 
the convergence. The 2-point statistics, such as the power 
spectra, therefore cannot contain all of the information in 
the map. Further cosmological information is contained 
in the higher order moments, however the number of de- 
grees of freedom rapidly becomes large making it difficult 
to present all of the information. For this reason we shall 
present only the 'collapsed' N-point functions, and since 
the non-Gaussianity is most clearly manifest in real (rather 




Fig. 4. — (top) Comparison of the convergence angular power spec- 
tra from the simulations at varying levels of grid resolution with the 
semi-analytic prediction computed using the method of Peacock & 
Dodds 119961 and Eq. 1171 . The spectra are averages from 10 real- 
izations of our high resolution simulation, with A 2 . = 1(1 + l)Ci/27T, 
for grids of 1024 , 2048 2 , and 4096 2 points, and a side length of 
300 7i — ^^Mpc. (bottom) The ratio of the spectra from our simula- 
tions with the semi-analytic prediction. 
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Fig. 5. — (top) The 3-dimensional matter power spectrum of 3 
boxes at selected redshifts for our high resolution simulation, com- 
pared against the s emi-an alytic prediction using the prescription of 
Peacock & Dodds 119961 . (bottom) The ratios of the simulated 
spectra with the semi-analytic prediction. 
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Table 1 

The skewness and kurtosis reported by several 
groups using simulations and analytical predictions 
(see text). All are for Sl ro = 0.3 = 1 — Qa and h = 0.7, 
and with 54 given in units of 10 4 . the values for v.w 

are for a semi-analytic interpolation between the 
predictions of hyper- extended perturbation theory 
(hept) and the perturbation theory result, while 
tj report values using the halo model approach. 
All other values are for simulations. 



than Fourier) space we shall present only real-space re- 
sults. In FigureHJwe show the skewness and kurtosis of the 
convergence on an angular scale 9. These are obtained by 
smoothing the k maps using the IDL command SMOOTH, 
which performs a real space boxcar average using square 
apertures with an integer number of pixels on a side. We 
then calculate the skewness 



Ss(6) 



Ml 



and the kurtosis 



S 4 (0) = 



(«2)-3(«g) s 



(18) 



(19) 



by averaging over the pixels in a given field and then av- 
eraging the results for ten realizations. Individual realiza- 
tions vary substantially, suggesting error bars of ~ 5% for 

53 and ~ 10% for 54, in line with the results of White & 
Hu (2000). Figure [7| shows the dependence on scale of the 
lower order moments for two cosmological models which 
differ only in the clustering strength, erg. 

An accurate measurement of the higher order moments, 
^3 and S4., from lensing maps may provide important 
constraints on the cosmological parameters f2 m and as 
(Bernardeau, van Waerbeke, & Mellier 1997; Jain & Sel- 
jak[1997). In perturbation theory, the variance of the con- 
vergence exhibits a strong dependence on both quantities, 
while the skewness and kurtosis are independent of the lat- 
ter. Although the independence of S3 and S4 from as is 
not expected to hold on scales of interest where the lens- 
ing signal can be measured, a joint measurement of these 
quantities should still be useful for breaking the degener- 
acy between J7 m and as if the actual dependence can be 
understood. This has motivated efforts to extend skewness 
and kurtosis predictions to smaller angular scales. 

We compare some of these calculations with a variety 
of numerical results to assess the degree of convergence 
between groups and methods. Table ^ presents S3 and 

54 reported by several groups, both from simulations and 
by analytic means. All are for cosmological models with 
Q m = 0.3 = 1 — f^A and h = 0.7. The moments from 
simulations are from JSW, HM (Hamana & Mellier 12*00 1(1 . 
WH (White & Hump, H2 (Hamana l2TMjl . and this 
paper (VW), while the last two rows are analytic predic- 
tions. The first of these, v.W (van Waerbeke 12000b) ). is 
for a semi-analytic fit to N-body simulations which tries 
to interpolate between the Hyper-Extended Perturbation 
Theory (HEPT) prediction (Scoccimarro & Frieman[1999 
Hui 1999) and the perturbation theory result. The second, 
TJ (Takada & Jain 2005]), uses the Halo model approach. 

We note that the skewness in the highly non-linear 
regime is predicted by HEPT to vary as a$ 0A . Lowest 
order perturbation theory predicts that S3 and S4 would 
be independent of erg. If we adjust all of the results to 
as = 0.9 the disagreements are lessened (for example us- 
ing HEPT S 3 (l') would scale to 144 for VWa and to 148 
for VWb). Keeping in mind that the simulation results 
have uncertainties of 10% or so, the agreement between 
groups is quite good, though more work needs to be done 
if we are to obtain precision cosmology from the higher 
order moments. 
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Fig. 6. — (bottom) The averaged skewness 53 and (top) kurtosis 
54 of the convergence for ten realizations of our high resolution 
simulation, as a function of smoothing scale 8. These are computed 
for two different values of inter-plane spacing Ax for grids of 4096 2 
(solid, dashed) and 2048 2 (dotted). As in Figure |7| the kurtosis is 
in units of 10 
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Fig. 7. — The ten realization average (bottom) variance, (re 2 ), 
(middle) skewness, S3, and (top) kurtosis, S4, of the convergence 
for our two high resolution simulations as a function of smoothing 
scale 8. For display purposes, £4 is in units of 10 4 , and (re 2 ) is in 
units of 10 -4 . 





#/deg 


0.05 


485 


0.10 


164 


0.15 


69.0 


0.20 


34.1 


0.25 


18.2 


0.30 


10.3 


0.40 


3.65 


0.50 


1.52 



Table 2 

The average number of peaks in convergence (with 

K > K m in) PER 1° X 1° FIELD. 



Finally it is interesting to look at the extrema of the lens- 
ing maps. We located the convergence peaks for our maps 
by selecting all pixels with k greater a selected threshold 
convergence K m i n , and then directly comparing the conver- 
gence of nearby pixels within a radius of 0.3 arcminutes. 
The largest value was then selected as the peak's center. 
Table|2shows the average number of peaks in a 1° x 1° field 
for selected values of K m i n . The average profile of these 
peaks is roughly consistent with that of a projected NFW 
profile (Bartelmann ll996l Navarro, Frenk, & White ITW7|) . 

5. NUMERICAL ISSUES 

The numerical convergence of the algorithm described 
in the previous section is limited by the size and resolu- 
tion of the N-body simulation and by the finite grids used 
to compute the distortion and deflection of the light rays. 
In the following section, we make use of the power spec- 
trum as a metric in order to test these limits under various 
conditions, and to evaluate the basic ray-tracing imple- 
mentation against the less approximate versions. We then 
comment on the numerical issues involving the skewness 
and kurtosis. 

5.1. N-body Simulations 

We use N-body simulations to calculate the evolution 
of the dark matter which governs the formation of large- 
scale structure. We use 3 simulations each of a ACDM 
model. The first simulation evolves 512 3 particles in a 
300 fr -1 Mpc bo x, usin g the TreePM-SPH code (see the ap- 
pendix of White[2002 ) but without any SPH particles. The 
cosmological model is chosen to provide a reasonable fit to 
a wide range of observations with fl m — 0.3, = 0.7, 
H = lOO/ikms^Mpc -1 with h = 0.7, n B h 2 = 0.02, 
n = 1 and erg = 1 (corresponding to Sh — 5.3 x 10~ 5 ). 
The simulation was started at z = 60 and evolved to the 
present with the full phase space distribution dumped ev- 
ery 100/i _1 Mpc from z ~ 2 to z = 0. The gravitational 
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Fig. 8. — As in the top panel of Figure^ except the convergence 
spectra from our high resol utio n simulation are now compared with 
those computed using Eq. I2UI . 
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force softening was of a spline form, with a "Plummer- 
equivalent" softening length of 20/i _1 kpc comoving, and 
the particle mass is 1.7 x 10 10 Ii~ 1 Mq. The 3-d mass power 
spectrum of this simulation is is given for selected redshifts 
in Figure |3J 

We use this simulation as our 'fiducial' model unless oth- 
erwise specified. To test the dependence on the amplitude 
of the power spectrum we also used a second simulation 
which is identical to the first except that the amplitude 
of clustering is reduced: eg = 0.8. We also make use of 
a smaller, lower resolution simulation for certain of the 
numerical tests described below. This final simulation is 
similar to the larger ones except that it uses fewer (128 3 ) 
particles and a smaller (128 /i _1 Mpc) box and it dumps the 
output more frequently (every 25 /i _1 Mpc). The softening 
length in this simulation is 36 /i _1 Mpc. 

5.2. Large scale resolution 

Our ability to resolve large scale features is limited by 
the field of view of our ray tracing, and therefore by the size 
of the box of the N-body simulation and the distance to 
the furthest lens plane. At z = 1, this limit is 7.5 degrees 
for our larger simulation and 3.2 degrees for our smaller 
one. However, most of our ray tracing runs were made at 
a field of view of 2 degrees. Only a few modes are avail- 
able to contribute to the power spectrum on large scales 
for any given realization, leading to substantial variation 
for different realizations. We therefore generate multiple 
realizations of ray tracing for an N-body simulation and 
average them to obtain estimates on large scales. For each 
realization, we project mass randomly along one of three 
different axes, and then select the origin of the projected 
plane at random. Depending on the distance from the ob- 
server, only a fraction of this projected plane is used in 
any given ray tracing realization. For example, a 2 degree 
run on our larger simulation uses less than 2% of a lens 
plane at the peak of the lensing kernel. Multiple nearly 
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Fig. 9. — As in Figure Isl except for our low resolution N-body 
simulation, which uses 128^ particles and has a box side length of 
128h _1 Mpc. The spectra are for grids of size 1024 2 , 2048 2 , and 
4096 2 . 



independent realizations are therefore possible using this 
method, and the results can be averaged to obtain lensing 
statistics. 

5.3. Small scale resolution 

The discrete nature of our ray tracing procedure intro- 
duces numerical smoothing when we interpolate mass den- 
sity values from the N-body simulation onto the Fourier 
grid. The N-body simulation is characterized by both a 
finite spatial resolution and by shot noise due to the finite 
particle number. To investigate the effect of this limita- 
tion on the angular resolution of the power spectrum, we 
model the finite resolution effect with Gaussian filters and 
thereby modify Eq. fTT| so that A^(fe) is replaced with 

Al(k) = (^Al(k)e-^ k2 + ^) e~^ (20) 

where A^(fc) is an effective 3-d mass power spectrum, 
fc 3 /27r 2 n is the analytic expression for the shot noise, n 
is the mean particle number density, and a n and o~ g are 
characteristic resolution limits of the N-body simulation 
and the Fourier grid, respectively. We did not find it nec- 
essary to introduce a term for shot-noise on each individual 
lens plane. 

We generate spectra using Eq. H20(l . normalizing the 
height to equal that of the simulation at i = 300. As 
shown in Fig. |H1 the spectra generated in this manner us- 
ing a g = 0.54Lbox/-^Vgrid (corresponding to a full width 
half max resolution of 1.3 grid cells) and a n set equal to 
zero are a good fit with those from our high resolution sim- 
ulation. This suggests that for this simulation the power 
spectrum resolution is limited by the grid and not by the 
N-body resolution. It is also consistent with our expecta- 
tion that the Fourier grid's resolution is roughly equal to 
its spacing. 

The limits imposed by the N-body simulation can be 
seen in the spectra of our low resolution simulation (Fig.^Jl, 
which are well approximated by those generated using 
Eq. ipD)l. with o- g — 0.54Lb ox /Agrid as before, and a n = 
50/i _1 kpc, or 5% of the mean inter-particle spacing n~3 . 
The three curves are for different values of N. The curve 
for N — 1024 is similar to those of Figure HO however the 
roll over is now caused by both the grid resolution and the 
N-body resolution. The curves at N = 2048 and N = 4096 
show substantially more power at high angular resolution. 
This is because the noise term, which is larger than the 
signal from this N-body simulation for £ ~ 10 4 and above, 
is no longer suppressed by the grid resolution. 

The value for a n is independent of grid spacing and is 
instead related to the resolution of the simulation. How- 
ever, it is not equal to the "Plummer equivalent" softening 
length, which is 36/i~ 1 kpc. To check this, we computed 
the 3-d mass power spectrum A^(fc) using the rebinning 
technique of Smith et al. (19999 for three simulations with 
identical initial conditions but varying particle number. 
We find that a n = 0.05 n~3 does roughly reproduce the 
roll off in power between the models, though it clearly 
doesn't get the whole form right because the models have 
more power at high-fc than this filter suggests. 
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We make use of our low resolution simulation to inves- 
tigate the convergence of the algorithm with respect to 
the number of lens planes used and to determine the ef- 
fect of decreasing the redshift/time interval between parti- 
cle dumps of the N-body simulation. Figure HU1 shows the 
power spectrum, averaged over 10 realizations, for different 
values of inter-plane spacing (Ax) and redshift/time inter- 
val (At). A change in Ax from 125/i _1 Mpc to 25/i~ 1 Mpc 
produced a change in the power spectrum roughly at the 
10% level for this simulation (for our high resolution sim- 
ulation this value was ~ 5%), however, additional de- 
creases in Ax did not significantly improve the resolu- 
tion. For a similar change in At, the spectra are within 
1% of each other. This is not surprising, since structures 
evolve on a scale roughly equal to the Hubble time, which 
is much longer than the time it takes for light to travel 
lOO/i^Mpc. 

We also wish to understand the numerical resolution 
of the weak lensing algorithm as it relates to the higher 
order moments. Figure [S] shows the effect on S3 and S4 
of changing the Fourier grid resolution and the inter-plane 
spacing A\. Increasing the resolution of the grid increases 
the small angle power in both the skewness and kurtosis, 
while decreasing A\ increases their power over the full 
range of resolution. However, this effect seems to be due 
entirely to the decrease in the convergence variance (k 2 ), 
while (k 3 ) and (k 4 ) are essentially unaffected. 

5.4. Discontinuities 



"angle" method described above. We randomly select the 
origin and orientation of the first box and then trace the 
rays through all the stacked boxes without any further ro- 
tations or origin shifts. The rays are then initialized in a 
field of view about a central line of sight that is angled with 
respect to the boxes to avoid repeating any structures, and 
the lens planes are made so that they are perpendicular 
this line of sight. We averaged the power spectra of ten 
runs created using this method and find that it is the same 
as for the standard "rotate and shift" method to better 
than 1%. 

5.5. 3- dimensional ray tracing 

One approximation inherent in the lens plane method 
is that since the light rays are all traveling in different 
directions they are only roughly perpendicular to the lens 
planes. The weak lensing formalism is based on Eq. QJ, 
da = — 2Vj_0 dXi which requires that we use the gradient 
perpendicular to the path of the light-rays. 

We used the 3-dimensional ray-tracing algorithm de- 
scribed previously to test the validity of this approxima- 
tion, and we find that it holds extremely well for two de- 
gree fields. We used three-dimensional Fourier grids with 
512 x 512 x 32 points, with the latter number applying to 
the direction of the line of sight. We find that convergence 
maps made with the two different methods agree to bet- 
ter than 0.1%, and measures such as the power spectra, 
skewness, and kurtosis are virtually indistinguishable. 



The process of rotation and origin shifting for each box 
necessarily leads to discontinuities in the 3-dimensional 
mass distributions, which may, in principle, introduce ar- 
tifacts during ray-tracing runs. To investigate this pos- 
sibility we traced through our large simulation using the 
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Fig. 10. — (top) The convergence power spectrum (averaged for 10 
realizations) of our low resolution simulation using different values 
of comoving distance Ax between lens planes and time intervals 
At between particle dumps of the N-body simulation. The units 
in both cases are h~ 1 Nlpc. (bottom) The ratio of the two spectra 
with Ax = 125 h~ 1 Mpc with respect to the one with Ax = At = 
25h~ 1 Mpc. 



6. PHYSICAL ISSUES 

6.1. The Born approximation near peaks in the 
convergence 

An approximation that is often useful in weak lensing 
analysis is the Born approximation, in which the lensing 




Ratio 

Fig. 11. — A histogram of the ratio k/kj of the convergence for 
300 large (re > 0.35) peaks in our maps. Here, k is the conver- 
gence at the center of a given peak using our full simulation, and 
Kg is the convergence of the same peak calculated using the Born 
approximation. 
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events for a given light ray are completely decoupled from 
one another. The distortion and deflection contributions 
for each event are computed at points along an undeflected 
path and as if for an undistorted ray. These are then 
summed to produce the total deflection and distortion. 
The Born approximation can be calculated by tracing rays 
along straight paths without deflections and approximat- 
ing A„ from Eq. Ijlljl as 

n-l 

A„ = I - ^ 9(Xm, Xn)U,n (21) 
m— 1 

The Born approximation has been shown before to be 
valid for whole convergence maps using measures such as 
the power spectrum. However, a breakdown of this ap- 
proximation, if one were to occur, would almost certainly 
take place where lensing effects are strongest, which could 
in principle be an issue of relevance to lensing by clusters 
of galaxies in the real universe. In Fig.^] we compare 300 
peaks with convergence values of at least 0.35 computed 
using our standard ray-tracing vs. the same peaks using 
the Born approximation. 

We find that the Born approximation is nearly always 
accurate to within 10% even for these large peaks, and the 
single biggest departure was 20%. In addition, the mean 
of the ratios is very nearly equal to one, so deviations ap- 
pear to be unbiased. We expect that the Born approxima- 
tion becomes increasingly accurate as k is decreased, but 
matching the 'Born peaks' with those of the full ray trac- 
ing (including deflections) becomes increasingly difficult so 
we have not been able to quantify this. 

6.2. Comparison to a universal magnification PDF 

In Figure |2 we compare the magnification probability 
distribution function of a typical run, P(n), with the ana- 
lytical prediction using the universal probability distribu- 
tion function P(?y) (Wang, Holz, & Munshi 2002), where 

P(M) = P(r?)/2| Kmin | . (22) 

The reduced convergence is defined as 77 = | (p— 1)/ 1 K m in | + 
1, and K m in is the "empty beam" convergence. 

The top panel shows a typical PDF from our simulations 
vs. the prediction from the model. The shape of the model 
PDF agrees with the simulation value to within ~ 5% near 
the PDF peak and ~ 20% in the high /x tail. However, we 
found that the model PDF was often shifted (in either 
direction) on the horizontal axis with respect to that of 
the simulation. The magnification /1 at the PDF peaks was 
typically within 0.5% but disagreed by as much as 1.5% 
for extreme cases, such as the one shown here. This is a 
substantial disagreement given the few percent amplitude 
of the lensing signal. Even so, the shape of the PDF is still 
well fit, as can be easily seen if the model PDF is shifted 
to lower n so that the curves overlap. 

6.3. Multiple lensing events 

During the ray-tracing runs, we tracked the number of 
lensing events above selected threshold magnitudes K m for 



each ray. A lensing event is defined here as occurring when 
a ray experiences a shear on a single lens plane that causes 
a change in the convergence at the observer in excess of 
that ff(Xm,X«)U m > n m in Eq. iJTTJ- In Tabled 
we present the results from ten simulations with a total of 
more than 4 x 10 7 lines of sight. It is quite uncommon for 
a given light-ray to experience multiple lensing events of a 
large magnitude, which suggests that a given large peak in 
the convergence maps is nearly always associated with a 
single localized massive object in the N-body simulation. 

7. CONCLUSIONS 

Weak gravitational lensing has become a powerful tool 
for observing the large-scale matter distribution in the 
universe, with rapid advances in observational capabili- 
ties hinting at even greater things to come. As observers 
achieve ever greater precision, it's important that the the- 
ory keep pace. 

In this paper, we have described a multi lens-plane algo- 
rithm for generating maps of weak lensing distortion from 
structure generated by N-body simulations paying partic- 
ular attention to the numerical convergence of the algo- 
rithm^). While our findings are substantially in line with 
previous works, we have noticed that biases can creep in if 
not enough lensing planes or time dumps are used. As con- 
cerns the range of validity of the generated maps, we find 
that the small-scale spatial resolution is well described by 
Eq. l|2"0"|l . With modern, parallel, codes the N-body simu- 
lation is not a limiting factor in lensing studies, suggesting 
that grids of models could be run relatively inexpensively. 

The basic multi lens-plane algorithm introduces two 
numerical artifacts to the calculation of distortion maps 
(other than finite resolution of the grids used in computa- 
tions) . One is the introduction of structure discontinuities 
at the box boundaries in the 3-dimensional matter distri- 
bution used to make the lens planes, and the other is that 
light-rays are not truly perpendicular to the lens planes. 
We have examined the effect of both of these artifacts using 
techniques specifically designed to avoid them, by tracing 
through the boxes at an angle in the first case, and us- 
ing a 3-dimensional simulation in the second, and we have 
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0.01 


40.8 


9.78 


1.61 


0.19 


0.05 


3.46 


0.06 


0.0005 


< 10- 


0.10 


0.68 


0.003 


< 10- 4 





0.15 


0.21 


0.0003 








0.20 


0.08 


0.0001 








0.25 


0.03 


< 10- 4 









Table 3 

The number of lensing events above a threshold 
convergence k m for 10 runs, given here as a 
percentage of the roughly 42 million total number 
of lines of sight. llght-rays are unlikely to 
experience large distortions more than once, 
suggesting that peaks in the convergence are 
almost always the result of a single lensing event. 
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shown that they do not introduce significant changes in 
the maps. 

We examined the applicability of the Born approxima- 
tion to weak lensing using measures of the statistical prop- 
erties of whole maps, such as the power spectrum and 
the skewness, and confirm that they are in generally good 
agreement with semi-analytic predictions. We also confirm 
that the Born approximation is valid even near peaks in 
the convergence, where the most non-linear lensing events 
in our simulations are likely to occur. 

Weak gravitational lensing is a tool that is providing 
insight into the nature of our universe at a rate that is 
expected to increase dramatically in coming years, and 
simulations will play a useful role. We have demonstrated 
that simulations of weak lensing based on a simple im- 
plementation of the multiple lens-plane algorithm gener- 
ate maps of lensing distortion that are a good approxima- 
tion of maps created using more computationally intensive 
methods. We have also shown that the numerical resolu- 
tion of these simulations is well approximated by a simple 
formula. This suggests that a relatively modest investment 
in computation can provide the highly accurate, simulated 
maps which are crucial to developing algorithms and un- 
derstanding data in this field. 

The simulations used here were performed on the IBM- 
SP2 at the National Energy Research Scientific Comput- 
ing Center. This research was supported by the NSF and 
NASA. 
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